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Abstract: We derive closed formulas for the condition number of a linear 
function of the total least squares solution. Given an over determined linear 
system Ax = b, we show that this condition number can be computed using 
the singular values and the right singular vectors of [A, b] and A. We also 
provide an upper bound that requires the computation of the largest and the 
smallest singular value of [A, b] and the smallest singular value of A. In numerical 
examples, we compare these values and the resulting forward error bounds with 
the error estimates given in |17j . 
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Une contribution au conditionnement du 
probleme de moindres carres totaux 

Resume : Nous derivons des formules exactes pour le conditionnement d'une 
fonction lineaire de la solution d'un probleme de moindres carres totaux. Etant 
donne un systeme lineaire sur-determine Ax = b, nous montrons que ce condi- 
tionnement peut etre calcule en utilisant les valeurs singulieres et les vecteurs 
singuliers a droite de [A, b] et A. Nous proposons aussi une borne superieure 
qui necessite le calcul de la plus grande et de la plus petite valeur singuliere 
de [A, b] ainsi que de la plus petite valeur singuliere de A. Dans des exemples 
numeriques, nous comparons ces conditionnements et les bornes d'erreur directe 
correspondantes avec les erreurs estimatives donnees dans |17| . 

Mots-cles : moindres carres totaux, conditionnement, perturbations norm- 
wise, modele d'erreur dans les variables 
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1 Introduction 

Given a matrix A £ l mx " (m > n) and an observation vector b £ M m , the 
standard over determined linear least squares (LS) problem consists in finding 
a vector x £ W 1 such that Ax is the best approximation of b. Such a problem 
can be formulated using what is referred to as the linear statistical model 

b = Ax + e, A £ R mxn , b £ R"\ rank(A) = n, 

where e is a vector of random errors having expected value E(e) = and 
variance-covariance V(e) = a 2 I. 

In the linear statistical model, random errors affect exclusively the observa- 
tion vector b while A is considered as known exactly. However it is often more 
realistic to consider that measurement errors might also affect A. This case is 
treated by the statistical model referred to as Errors-In- Variables model (see 
e.g [T7J p. 230] and [5J p. 176]), where we have the relation 

(A + E)x = b + e. 

In general it is assumed in this model that the rows of [E, e] axe independently 
and identically distributed with common zero mean vector and common covari- 
ance matrix. The corresponding linear algebra problem, discussed originally 
in |12j . is called the Total Least Squares (TLS) problem and can be expressed 
as: 

mm\\(E,e)\\ F , {A + E)x = b + e 1 (1) 

where ||-||^ denotes the Frobenius matrix norm. As mentioned in [T71 p. 238], 
the TLS method enables us to obtain a more accurate solution when entries of 
A are perturbed under certain conditions. 

In error analysis, condition numbers are considered as fundamental tools 
since they measure the effect on the solution of small changes in the data. In 
particular the conditioning of the least squares problem was extensively studied 
in the numerical linear algebra literature (see e.g [FJ [El EE [TO1 HFJ [TBI HFJ HI 
|2"2"]). The more general case of the conditioning of a linear function of an LS 
solution was studied in [5] and [U when perturbations on data are measured 
respectively normwise and componentwise (note that the componentwise and 
normwise condition numbers for LS problems were also treated in [9] but without 
the generalization to a linear function of the solution). Moreover we can find 
in [3J algorithms using the software libraries LAPACK [T] and ScaLAPACK [5J 
as well as physical applications. 

The notion of Total Least Squares was initially defined in the seminal pa- 
per [12] that was the first to propose a numerically stable algorithm. Then var- 
ious aspects of the TLS problem were developed in the comprehensive book [T7J 
including a large survey of theoretical bases, computational methods and ap- 
plications but also sensitivity analysis with for instance upper bounds for the 
TLS perturbation. The so-called Scaled Total Least Squares (STLS) prob- 
lem (mhi£ je || [E, e)||jr, (A + E)xj = 76 + 6, for a given scaling parameter 7) 
was formulated in [20] in which were addressed the difficulties coming from non 
existence of TLS solution. In a recent paper [23] , we can find sharp estimates of 
the normwise, mixed and componentwise condition numbers of the Scaled Total 
Least Squares (STLS) problem. 
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Here we are concerned with the TLS problem, which is a special case of 
the STLS problem, and we will consider perturbations on data (A, b) that are 
measured normwise using a product norm. Contrary to [23] . we will consider 
the general case of the conditioning of L T x, linear function of the TLS solution 
for which we will derive an exact formula. The common situations correspond 
to the special cases where L is the identity matrix (condition number of the 
TLS solution) or a canonical vector (condition number of one solution compo- 
nent). The conditioning of a nonlinear function of a TLS solution can also be 
obtained by replacing, in the condition number expression, the quantity L T by 
the Jacobian matrix at the solution. 

We notice that the normwise condition number expression proposed in |23j 
is based on the evaluation of the norm of a matrix expressed as a Kronecker 
product resulting in large matrices which may be, as pointed out by the authors, 
impractical to compute, especially for large size problems. We propose here a 
computable expression for the resulting condition number (exact formula and 
upper bound) using data that could be already available from the TLS solution 
process, namely by-products of the SVD decomposition of A and [A, b}. We also 
make use of the adjoint operator which enables us to work on a space of lower 
dimension and we propose a practical algorithm based on the power method. 



2 Definitions and notations 
2.1 The total least squares problem 

Let A e R mx " and b G R™, with m > n. Following [T7], we consider the 
two singular value decompositions of A, and [A, b] : A — U'H'V T and [A, b] = 
UT,V T . We also set £ = diag(eri, . . . , er„ +1 ), £' = diag(a^, . . . , a' n ), where the 
singular values are in nonincreasing order, and define A,; = of, and A J = a?. 
From [5, p. 178], we have the interlacing property 

01 > 0i > cr 2 > • •■ > a n > a' n > <r n+1 . (2) 

We consider the total least squares problem expressed in Equation |T]) and 
we assume in this text that the genericity condition a' n > ctVi+i holds (for 
more information about the "nongeneric" problem see e.g |17l I20| ). From [TT1 
Theorems 2.6 and 2.7], it follows that the TLS solution x exists, is unique, and 
satisfies 

x= (A T A-X n+1 I n y 1 A T b. (3) 



In addition, 



-1 



is an eigenvector of [A, b] T [A, b] associated with the simple 



eigenvalue A„ + i, i.e a' n > a n +i guarantees that A n+ i is not a semi-simple eigen- 
value of [A, b] T [A, b}. As for linear least squares problems, we define the total 
least squares residual r = b — Ax, which enables us to write 



A 



ra+l 



1 + X 1 X 



■1] 



' A T A 


A T b 




X 


b T A 


b T b 




-1 



1 + X 1 X 



(4) 



As mentioned in [T71 p. 35], the TLS solution is obtained by scaling the last 
right singular vector of [^4, b] until its last, component is —1 and, if Wi jrl +i 
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denotes the ith component of v n +i, we have 



X= [Vl, n +1, ■ ■ . ,«n,"+l] T - ( 5 ) 

V n +l,n+l 

The TLS method involves an SVD computation and the computational cost is 
higher than that of a classical LS problem (about 2mn 2 + 12n 3 as mentioned 
in [131 p. 598], to be compared with the approximately 2mn 2 flops required 
for LS solved via Householder QR factorization). However, there exist faster 
methods referred to as "partial SVD" (PSVD) that calculate only the last right 
singular vector or a basis of the right singular subspace associated with the 
smallest singular values of [A, b] (see [TTJ p. 97]). 

2.2 Condition number of the TLS problem 

To measure the perturbations on data A and b, we consider the product norm 
defined on R m x " x R m by 1 1 ( A, b) \ \ F = ^\\A\\ 2 F + \\bf 2 and we take the Euclidean 
norm ||x|| 2 for the solution space M. n . In the following, the n x n identity matrix 
is denoted by /„. 

Let L be a given n x k matrix, with k < n. We suppose here that L is not 
perturbed numerically and we consider the mapping 

g : W' nx n x W n — > R k 

(A,b) i v g(A,b)=L T x = L T (A T A-\ n+1 I n )- 1 A T b, 

Since A„+i is simple, g is a Frechet-differentiable function of A and b, and the 
genericity assumption ensures that the matrix (A T A— A„ + i/„) _1 is also Frechet- 
differentiable in a neighborhood of (A, b). As a result, g is Frechet-differentiable 
in a neighborhood of (A, b). 

The approach that we follow here is based on the work by [TTJ [2T] where the 
mathematical difficulty of a problem is measured by the norm of the Frechet 
derivative of the problem solution expressed as a function of data. This measure 
is an attainable bound at first order, and may therefore be approximate when 
large perturbations are considered. 

Using the definition given in [111 I21| , we can express the condition number 
of L T x, linear function of the TLS solution as 

KIT A h\ ||ff'(Afr)-(AA,Afe)|| 2 

K(L,A,b)= max ,. . - - — — — -. (6) 

K(L,A,b) is sometimes called the absolute condition number of L T x as op- 
posed to the relative condition number of L T x and defined, when L T x is nonzero 

by 

if W (i, A, b) = K(L, A, 6)|| (A, b)\\ F / \\L T x\\ 2 . (7) 

In the remainder, the quantity K(L,A, b) will be simply referred to as the 
TLS condition number, even though the proper conditioning of the TLS solution 
corresponds to the special case when L is the identity matrix. 
In the expression g'(A, 6).(A^4, A6), the "." operator denotes that we apply the 
linear function g'(A, b) to the variable (AA,Ab). We will use this notation 
throughout this paper to designate the image of a vector or a matrix by a linear 
function. 
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Remark 1 The case where g(A,b) = h(x), with h being a differentiable non- 
linear function mapping M. n to R fc is also covered because we have 

g'(A, b).(AA, Ab) = h'(x).(x'(A, b).(AA, Ab)), 

and L T would correspond to the Jacobian matrix h'(x). The nonlinear function 
h can be for instance the Euclidean norm of part of the solution (e.g in the 
computation of Fourier coefficients when we are interested in the quantity of 
signal in a given frequency band). 



3 Explicit formula for the TLS condition number 

3.1 Frechet derivative 

In this section, we compute the Frechet derivative of g under the genericity 
assumption, which enables us to obtain an explicit formula for the TLS condition 
number in Proposition [5] 

Proposition 1 Under the genericity assumption, g is Frechet differentiable in 
a neighborhood of (A, b). Setting B\ = A T A — X n+ il n , the Frechet derivative of 
g at (A, b) is expressed by 

g'(A, b) : R mxn x W n — > R k 

(AA, Ab) i — > L T B^ (a t + ^^) (Ab-AAx)+ (8) 

L T B- x x AA T r. 

Proof: The result is obtained from the chain rule. Since A n +i , expressed in Equa- 
tion (H]), is a simple eigenvalue of [A, 6] T [A, b] with corresponding unit eigenvec- 
tor , 1 \x T — ll , An+i is differentiable in a neighborhood of (A, b) and 

\l \ J rX T x 

then we have 

y n+1 (Ab).(AA,Ab) = — V[/ -1] 



1 


1 - 


h X T X 




2 


1 - 


h x T x 




2 


1 - 


h X T X 




2 


1 - 


h X T X 



AA T A + A T AA AA T b + A T Ab 
b T AA + Ab T A Ab T b + b T Ab 





X 




-1 



yielding 



[x T AA T Ax ~ x T AA T b - x T A T Ab + b T Ab) 
(-x T AA T r + (b T - x T A T )Ab) 
(-r T AAx + r T Ab) , 



/a 7\/a a A7\ 2r T (Ab — A Ax) 

X' n+1 (A, b).(AA, Ab = y — >-. 9 

1 + x 1 x 



Applying the chain rule to B x 1 , we obtain 

T A , at , 



(B^)'(A, b).(AA, Ab) = -B- 1 (AA T A + A T AA - X' n+1 (A, b).(AA, Ab)I n ) B^ 1 

<Ta , aT AA 2r T (Ab-AA x) 

1 + X T X 



-B- 1 ( AA T A + A T AA - - I n ) B y 
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The chain rule now applied to g (A, b) leads to 

g'(A, b).(AA, Ab) = -L T B^ (AA T A + A T AA - X' n+1 (A, b).(AA, Ab)I n ) B^ 1 A T b + L T B^ 1 AA T b + L T . 

= -L T B^ (AA T A + A T AA - X' n+1 {A, b).{AA, Ab)I n ) x + L T B^ 1 (AA T b + A T Ab) 
2xr T 



A 1 



1 + x T x 



(Ab — AAx) + L 1 B7 L AA I r, 



which gives the result. 



□ 



We now introduce the vec operation that stacks all the columns of a matrix 



into a long vector: for A = [<n,...,a n ] <E 



vec(,4) = [a(,...,a^] T g 
E" 1 " xl . Let P e l mnXm " denote the permutation matrix that represents the 
matrix transpose by vec(B T ) = Pvec(B). We remind also that vec(AXB) = 
(B T Cg) A)vec(X), where <£) denotes the Kronecker product of two matrices |141 
p. 21]. 

Let us now express the matrix representing g'(A, 6), denoted by M g >. Since 
g'(A, b).(AA, Ab) e R k , we have g'(A, b).(AA, Ab) = vec(g'(A, b).(AA, Ab)) and 



setting in addition D\ 
g'(A,b).(AA,Ab) = 



A T 



1+X T X 



pkxm 



we obtain from ([5]) 



vec (D x (Ab - AAx) + L T B^AA T r) 

(-x T <g> D x ) vec(AA) + (r T <g> (L T B^ 1 )) vec(AA T ) + D x Ab 



= [-x T ® Dx 



(L T B^))P, Dx] 



vec(AA) 
Ab 



Then we get 

Mg, = [-X 



Dx+(r 
But we have \\(AA,Ab)\\ F 



^ (L T B^)) P, D x ] 6 R kx ( nm + m \ 

and then, from Proposition [T] and 



vec(AA) 
Ab 

using the definition of K(L,A,b) given in Expression ©, we get the following 
proposition that expresses the TLS condition number in terms of the norm of a 
matrix. 

Proposition 2 The condition number ofg(A,b) is given by 

K(L,A,b) = \\M g ,\\ 2 , 

where 



[-x T ®D X + (r T ® (L T Bl 1 )) P, D x ] & 



pkx (nm+m) 



3.2 Adjoint operator and algorithm 

Computing K(L, A, b) reduces to computing the spectral norm of the k x (nm + 
m) matrix M g '. For large values of n or m, it is not possible to build explic- 
itly the generally dense matrix M g >. Iterative techniques based on the power 
method [TH p. 289] or on the Lanczos method [12] are better suited. These 
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algorithms involve however the computation of the product of -M^, by a vector 
y € K. fe . We describe now how to perform this operation. 

Using successively the fact that B X T = P^ 1 , (A® B) T = A T B T ', 
vec(AXB) = (B T ® v4)vec(X) and P T = p- 1 we have 



M T g ,y 



-x®Dl + P T (r ® (B X T L)) 



y 

-1 ; 



-(a ® P^)vec(y) + P T (r ® (P A vec(y) 



£>A> 



P- 1 (Pvec (-Dlyx T ) + vec (B x 1 Lyr T )) 
Dlv 



P- 1 (vec ((-Dlyx T ) T ) + vec (B x 1 Lyr T )) 

|T 
A 



P _1 vec {-xy T D x + B x Lyr 7 
Dlv 



and since for any matrix P we have P 1 vec(B) = vec(B T ), we get 



vec (-P^ya; T + ry T L T B x n 



(10) 



This leads us to the following proposition. 

Proposition 3 TTie adjoint operator of g' '(A,b) using the scalar products 
tra,ce(Aj A 2 ) + bjb 2 and y\yi respectively on K mxn x K m and K fe is 



In addition, if k = 1 we /iGroe 



p m x n w Tip m 



(-D T x yx T + ry T L T B x \D T x y) 



(11) 



K{L,A,b) = J\\-Dlx T + rL T B 



A IIF 



I^aII* 



(12) 



Proof: Let us denote by < (A\, bx), (A2, 62) > the scalar product trace^f.^)- 

x 



6f 62 on jR mXTl x M m . We have for any y G M fe 
!/ r ( 5 '(A6).(A4,A&)) = y T M. g i 



vec(AA) 
A& 



vec(Avl) 
Afc 



= vec (-Dlyx T + ry T L T P^ 1 ) vec(AA) + (Dj y) T Ab. 

Using now the fact that, for matrices A\ and A2 of identical sizes, 
vec(Ai) T vec(j42) = trace (A^^), we get 



y T (g'(A, b).(AA, Aft)) = trace {(-D T x yx T + ry L B x ) AA) + [D T x y) T Ab 
= < {-D T x yx T + ry T L T B x \Dly) , (AA, Ab) > 
= < g'*(A, b).y, (AA, Ab) >, 
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which concludes the first part of the proof. 
For the second part, we use 



vec (—D x yx T + ry T L T B x 



K(L, A, b) = \\M g >\\ 2 = ||AV|| 2 =max 



vec(~D^x T + rL T B~ u 



vec(£>D, 



Since k = 1, we have y € R, and iT(Z/,A, &) = 
and the result follows from the relation vec(.<4i) T vec(.<4i) = traccA^Ai = 

WMl- 

□ 

Remark 2 The special case k = 1 recovers the situation where we compute the 
conditioning of the ith solution component. In that case L is the ith canonical 

J^ 1 is the ith row of B x 



vector of R" and, in Equation (fT2l . L T B X 1 is the ith row of B A 1 and £>a is the 
ith row of B7 1 (A T 



Using ((HJl and (|TT|) . we can now write in Algorithm 1 the iteration of the 
power method ( |16l p. 289]) to compute the TLS condition number K(L, A, b). 



Algorithm 1 : Condition number of TLS problem 
Select initial vector y £ M. k 
for p=l,2,... 

(A p ,b p ) = {-D T x yx T + ry T L T B x \D T x y) 

v = \\(A p ,bp)\\ F 

(A p ,b p ) <- (i • Ap,± -b p ) 

y = I?B? (A T + ^) (b p - A p x) + iFB^r 

end 

K(L, A,b) = y/v 



The quantity v computed by Algorithm 1 is the largest eigenvalue of Mg/M.^ 
Since K(L,A,b) = ||A4 g '|| 2 then the condition number K(L,A,b) is also the 
largest singular value of M g > i.e ^Jv. As mentioned in |131 p. 331], the al- 
gorithm will converge if the initial y has a component in the direction of the 
corresponding dominant eigenvector of Ai g 'A4^,. When there is an estimate 
of this dominant eigenvector, the initial y can be set to this estimate but in 
many implementations, y is initialized as a random vector. The algorithm is 
terminated by a "sufficiently" large number of iterations or by evaluating the 
difference between two successive values of v and comparing it to a tolerance 
given by the user. 
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3.3 Closed formula 

Using the adjoint formulas obtained in Section [3.21 we now get a closed formula 
for the total least squares conditioning. 

Theorem 1 We consider the total least squares problem and assume that the 
genericity assumption holds. Setting B x = A T A — A„+i/„, then the condition 
number of L T x, linear function of the TLS solution, is expressed by 

K(L,A,b) = \\C\\l , 
where C is the k x k symmetric matrix 

C = (1 + \\xf 2 )L T B^ (a t A + A„ +1 (7„ - -^f^)) B?L. 

V 1 + IMI 2 / 



Proof: We have K(L,A,b) 2 — ||A^p|| 2 = max|| y || 2=1 . If y is a unit 

vector in R fe , then using Equation ([TU]) we obtain 

\\M T g ,yf 2 = \\vec(-Dlyx T + ry T L T B x 1 )\\l + \\Dly\\l 
= \\-Dlyx T + ry T L T B- 1 \\ 2 F +\\Dly\\l 

= \\Dlyx T \\ 2 F + \\ry T L T B^\\ 2 F - 2tv & ce{xy T D x ry T L T B x x ) + \\D 
For all vectors u and v, we have ||ww T |ji? = ||w|| 2 |Mj 2 - Moreover we have 
trace ({xy T D x r){y T L T B^ 1 )) = trace ({y T L T B^){xy T D x r)) = y T L T B x 1 xr T D 
Thus 

\\M T g ,y\\l = \\x\\l \\Dly\\l + \\rf 2 \\B^Lyf 2 - 2 y T L T B x ' xr T D T x y + \\D T x y\\\ 
= (l + x T x)y T D x Dly+\\r\\ly T L T B- x 2 Ly-2y T L T B^xr T D T x y 
= y T ((1 + x T x)D x D T x + \\r\\l L T B X 2 L - 2 I^B^x^D^ y, 

i.e \\M%\\1 = \\C\\ 2 with 

C = (1 + x T x)D x Dl + \\r\\ 2 2 L T B X 2 L - 2 L T B x x xr T D T x . (13) 
Replacing D x by L T B x x (a t + 2 ^ x ) , Equation (T3]) simplifies to 

C = L T B X 1 ({l + x T x)A T A+\\r\\ll n + 2A T rx T ^B x 1 L. (14) 

But A T rx T = A T (b — Ax)x T = A T bx T — A T Axx T and, since from Equa- 
tion ^ we have A T b = B x x, we get A T rx T = B x xx T - A T Axx T = (A T A - 
Ki+iIn)xx T — A T Axx T = — \ n +ixx T . From Equation ^ we also have 1/ n ; 
A n +i(l + x T x) and thus Equation p^|) becomes 



C = L T Bl x {{l + x T x)A T A + X n+1 {l + x T x)I n ~ 2\ n+1 xx T ) B^L 



Y^a) ) B x l1 



□ 
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4 TLS condition number and SVD 
4.1 Closed formula and upper bound 

Computing K(L,A,b) using Theorem [T] requires the explicit formation of the 
normal equations matrix A T A which is a source of rounding errors and also 
generates an extra computational cost of about ran 2 flops. In practice the TLS 
solution is obtained by Equation ([5]) and involves an SVD computation. In the 
following theorem, we propose a formula for K(L,A, b) that can be computed 
with quantities that may be already available from the solution process. In the 
following ni i (resp. 0i iTl ) denotes the zero column (resp. row) vector of length 

71. 

Theorem 2 Let V and V' be the matrices whose columns are the right sin- 
gular vectors of respectively [A, b] and A associated with the singular values 
(o"i, . . . , ayi+i) and (a[, . . . , o' n ). Then the condition number of L T x, linear func- 
tion of the TLS solution is expressed by 



K(L,A,b) = (l 



|2\i 



L T V'D'[v' T , O n ,i]V[D, nA ] 



where 



D> = diag ((a; 2 - <7 2 +1 )-\ • ■ • , (<# - o^+i) -1 ) and D = diag ((a 2 + <J 2 n+1 )i , . . . , (a 2 + <r 2 n+1 )^ 
When L is the identity matrix, then the condition number reduces to 



K(L,A,b) = (l + \\x\\l)i D'[V' T , O nA ]V[D, nA ]- 



Proof: From [A, b] = UZV T , we have [A, b] T [A, b] = VY?V T = YTiZl °f*W 



and 



n+1 



n+1 



[A,b] T [A,b} + X n+1 L n+1 = ^ Q-fvivf + A„+i ^ VjvJ 

i=l i=l 
n+1 

i=l 
n 

= + ^n+lM^f + 2A n+l w ™+l u n+l> 



leading to 



i=l 



[A, 6] T L4, 6] + An+i/n+i - 2A„+iU„+i^ +1 = J^(cr 2 + al +1 )v t vf (15) 

and, since u„+i is a unit 



1 



From Equation ([5]), we have v n +i = — iVi+i,n+i 

vector, w 2 +1 „ +1 = x+iLlI" ' Then Equation (|15[) can be expressed in matrix 
notation as 



' A T A 


A T b 


+ A n +1 


In On,l 


2A n +i 


r t i 

XX —x 


b T A 


b T b 


0l,„ 1 _ 


i+h! 


-x T 1 



X>i+<^+i) 



1=1 



(16) 
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The quantity A T A + X n+ i(I n — corresponds to the left-hand side of 

Equation (|16p in which the last row and the last column have been removed. 
Thus it can also be written 



2xx r I v-^ 

A T A + A n+ i(/„ - — s ) = [I n , n ,i] }Jct 2 + f7 «+i) 7 

1+INIa \<=i 
and the matrix C from Theorem [T] can be expressed 





'In' 


) 


Pi, n . 



C = (1 + |M|*)L T [S- 1 , On.1] X)(of + 



i=l 



) 






0i,„_ 



(17) 



Moreover from A = U'YIV" 1 ', we have A T ,4 = y'£'V T = £" =1 a '? v 'i^ and 



A A — X n+ \I n 



n 



E'2 / 'T 2 / 'j 



i=l 



T 



Hence S^ 1 = F'^W'" 1 = V"£>'V' T and [B^ 1 , 0„,i] = K'D' [v' T , n ,i] . 

•2 I ^2 \„, ,,T 



We also have £ j=1 0f + fn+iM^ = V 



£) 

0i,„ 



[D, O n ,l]V q 



Then, by replacing in Equation (|17|) . we obtain C = (1 + j|a;|| 2 )y^ T with 
y = L T V'D' [V' T , n ,i] V [D, n ,i] T - As a result, using Theoremffl 



if(i,A&) 2 = l|C|| 2 = (i + INl2) 



vv 1 



(1 + NI2) 



r 



When L = /„, we use the fact that V is an orthogonal matrix and can be 

~ 2 

removed from the expression of V 

2 

□ 

In many applications, an upper bound would be sufficient to give an estimate 
of the conditioning of the TLS solution. The following corollary gives an upper 
bound for K(L, A,b). 

Corollary 1 The condition number of L T x, linear function of the TLS solution 
is bounded by 



K(L,A,b) = (l + \\x\\l)*\\L\\. 



{vl+o-l +l )- 

Wn -°"n+l) 



Proof: This result comes from the inequality ||AB|| 2 < ||j4|| 2 ||S|| 2 , followed by 
\\D'\\ 2 =max l ( < r.i 2 -cr r 2 l+1 )- 1 = (^ 2 -cr 2 +1 ) _1 and ||L>|| 2 = max^cr 2 + ct 2 +1 ) = 
K+^+i)- 



□ 
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4.2 Numerical examples 

In the following examples we study the condition number of x i.e L is here the 
identity matrix /„. Then, to simplify the notations, we removed the variable L 
from the expressions and the condition number of x will be denoted by K(A, b) 
and its upper bound by K(A, b). All the experiments were performed with 
MATLAB 7.6.0 using a machine precision 2.22 • 10" 16 . 



4.2.1 First example 

In the first example we consider the TLS problem Ax ~ b where [A, b] is defined 

by 

[A,b]=Y(® ^jZ T e R" IX (» +1 ), Y = I m - 2yy T , Z = J n+1 - 2zz T , 

where y € R m and z 6 R" +1 are random unit vectors, D = diag(n,n — 
1, • • • ,1,1 — e p ) for a given parameter e p . The quantity a' n — n +i measures the 
distance of our problem to nongenericity and, due to Equation we have in 
exact arithmetic 

c„ - O'n+l < On — &n+l = e p . 

Then by varying e p , we can generate different TLS problems and by considering 
small values of e p , it is possible to study the behavior of the TLS condition 
number in the context of close-to-nongeneric problems. The TLS solution x is 
computed using an SVD of [A, b] and Equation ([S]) . 

In Table [TJ we compare the exact condition number K {A, b) given in The- 
orem^ the upper bound K(A, b) given in Corollary [TJ and the upper bound 
obtained from |17l p. 212] and expressed by 



k(A, b) 



9cti Ibll 



O'n+l 



+ 1/ ll & ll 2 -° 



n+1 



We also report the condition number computed by Algorithm 1, denoted by 
K p {A, b), and the corresponding number of power iterations (the algorithm ter- 
minates when the difference between two successive values is lower than 10 -8 ). 
When a' n — tr rl +i decreases, the TLS problem becomes worse conditioned and 
there is a factor O(10) between the exact condition number K (A, b) and its 
upper bound K(A, b). We also observe that K(A, b) is an estimate of better 
order of magnitude than k(A, b) and that, for small values of a' n — o n +i, k{A, b) 
is much less reliable. K p (A, b) is always equal or very close to K(A, b). 

4.2.2 Second example 

Let us now consider the following example from |17l p. 42] also used in |23j 
where 



.4 



/ m- 1 -1 

-1 m- 1 



-1 -1 
-1 -1 

V -i -i 



-1 \ 

-1 



m — 1 
-1 

-1 / 



ix(m- 2) t 



/ "I \ 

-1 



-1 
TO — 1 

V -i J 
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Table 1: TLS conditioning for several values of a' n — <r n +i- 



a n - a n+l 


K(A,b) 


K(A,b) 


k(A,V) 


K p (A,b) 


#iter 


9.99976032- lO" 1 


1.18 -10 u 


2.36- 10 1 


1.29- 10 2 


1.18- 10 u 


11 


9.99952397- 10~ b 


8.36 • 10 3 


1.18- 10 5 


1.31 • 10 1U 


8.36 ■ 10 3 


6 


9.99952365 • 10~ y 


8.36 • 10 Y 


1.18 • 10 y 


1.31 • 10 1 * 


8.36 • 10 Y 


4 


9.99644811 • 10~ 13 


8.36- 10 11 


1.18- 10 13 


1.31 • 10 2B 


8.32- 10 11 


5 



The exact solution of the TLS problem Ax sa b can be computed analytically |17[ 
p. 42] and is equal to x = — (1, • • • , 1) T . We consider a random perturbation 
(AA, Ab) of small norm || (AA, Ab)\\p = 10~ 10 and we denote by x the computed 
solution of the perturbed system (A + AA)x «i + Ab. 

In Table [2 we report for several values of m the relative condition number 
as defined in Q and we compare the computed relative forward error ^^jj^ 2 
with the forward error bounds that can be expected from the computation of 
K^ rel \A, b) and its upper bounds K^ rel \A, b) and ^^{A, b). Since the condi- 
tion number corresponds to the worst case in error amplification at first order, 
these quantities are, as observed in Table [21 always larger than the computed 
forward error (there is approximately a factor 10 2 between those quantities). We 
also observe that, in this example, K^ rel \A, b) and K^ rel \A, b) produce forward 
error estimates that are of same order of magnitude. 



Table 2: Forward error and upper bounds for a perturbed TLS problem. 



m 


K^ rel \A, b) 


||x-x|| 2 


fr(ret)/ A u\ ll(AA,A6)||p 


ff(rel)/ A u\ ll(AA,A6)|| F 


r (rel)( A M ll(AA,A6)|| 


\\x\\r, 


[A ^> \\{AM\W 


( A ' b) \\(A,b)\\ F 


^ b > \\{A,b)\\ F 


50 


5.05 • 10 1 


2.45 ■ 10~ 13 


2.21 • 10~ n 


1.55- 10- 1U 


6.72- lO" 10 


100 


1.01 • 10 2 


1.08 • 10~ 13 


1.55- lO^ 11 


1.54- 10~ iU 


4.26- lO" 1 " 


500 


5.01 • 10 2 


8.79- 10- 14 


6.85- 10~ i2 


1.53- 10- 1U 


1.66- lO" 10 


1000 


1.00- 10 3 


4.33- 10~ 14 


4.84- 10~ 12 


1.53- lO" 10 


1.13- lO" 10 



5 Conclusion 

We proposed sensitivity analysis tools for the total least squares problem when 
the genericity condition is satisfied. We provided closed formulas for the con- 
dition number of a linear function of the TLS solution when the perturbations 
of data are measured normwise. We also described an algorithm based on an 
adjoint formula and we expressed this condition number and an upper bound 
of it in terms of the SVDs of [A, b] and A. We illustrated the use for these 
quantities in two numerical examples. 
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